Antje Hofmann: Learning Syntactic Relations with the Mole Task (2)

RePsychLing in SMLP2022

1 Background

This version removes the outlier subject “p10” found in the version 1.

Children’s (age: 4-8 years)reaction times in a task teaching them syntactic relations.

1.1 Overview

• Original analysis is by Antje Hofmann.
• MixedModels.jl version
• Addition of new chunks illustrate
  • selection of parsimonious LMM using random-effects PCA
  • plotting conditional means
  • illustration of borrowing strength

2 Readme

2.1 Variables

• Subj: Participant ID (renamed from ID; random factor)
• Item: Word ID (random factor)
• age: 4 - 8 years
• Block (within-Subj/within-Item):
  • 1st Learning
  • 2nd Learning
  • Disruption
  • Recovery
• Target(renamend fom targetness)
  • non-syllable target
  • syllable target
• rt: response time

3 Setup

3.1 Packages

First attach the MixedModels.jl package and other packages for plotting. The CairoMakie.jl package allows the Makie graphics system [@Danisch2021] to generate high quality static images. Activate that package with the SVG (Scalable Vector Graphics) backend.

Code

using AlgebraOfGraphics
using Arrow
using CairoMakie       # graphics back-end
using CategoricalArrays
using Chain
using DataFrames
using DataFrameMacros  # simplified dplyr-like data wrangling
using MixedModels
using MixedModelsMakie # diagnostic plots
using ProgressMeter
using Random           # random number generators
using RCall            # call R from Julia
using StatsModels

using AlgebraOfGraphics: boxplot
using AlgebraOfGraphics: density

using MixedModelsMakie: qqnorm
using MixedModelsMakie: ridgeplot
using MixedModelsMakie: scatter
using MixedModelsMakie: caterpillar

ProgressMeter.ijulia_behavior(:clear);
CairoMakie.activate!(; type="svg");

• The data are available as an arrow file.
• Most of preprocessing was done with R in RStudio (see Hofmann_Maulwurf.Rmd).
• Order of factor levels should be checked.

Code

dat = DataFrame(Arrow.Table("./data/Hofmann_Maulwurf_rt.arrow"))
transform!(dat,
    :Target => categorical => :Target,
    :Block  => categorical => :Block,
    :age => (x -> x .- 6) => :a1, # center age at six years
    :rt => (x -> log.(x)) => :lrt)
describe(dat)
dat = filter(row -> row.Subj != "p10", dat)

6,257 rows × 8 columns (omitted printing of 1 columns)

	Subj	Item	age	Block	Target	rt	a1
	String?	String?	Float64?	Cat…?	Cat…?	Float64?	Float64
1	p01	02SG.ogg	7.73	1st Learning	Non-target syllable	3285.0	1.73
2	p01	21SG.ogg	7.73	1st Learning	Non-target syllable	2873.0	1.73
3	p01	20SG.ogg	7.73	1st Learning	Non-target syllable	2169.0	1.73
4	p01	13SG.ogg	7.73	1st Learning	Non-target syllable	3133.0	1.73
5	p01	32SG.ogg	7.73	1st Learning	Non-target syllable	5793.0	1.73
6	p01	05SG.ogg	7.73	1st Learning	Non-target syllable	2064.0	1.73
7	p01	27SG.ogg	7.73	1st Learning	Non-target syllable	4546.0	1.73
8	p01	18SG.ogg	7.73	1st Learning	Non-target syllable	2949.0	1.73
9	p01	23SG.ogg	7.73	1st Learning	Non-target syllable	2957.0	1.73
10	p01	07SG.ogg	7.73	1st Learning	Non-target syllable	2508.0	1.73
11	p01	29SG.ogg	7.73	1st Learning	Non-target syllable	2232.0	1.73
12	p01	19SG.ogg	7.73	1st Learning	Non-target syllable	2517.0	1.73
13	p01	14SG.ogg	7.73	1st Learning	Non-target syllable	2234.0	1.73
14	p01	16SG.ogg	7.73	1st Learning	Non-target syllable	4272.0	1.73
15	p01	15SG.ogg	7.73	1st Learning	Non-target syllable	3389.0	1.73
16	p01	09SG.ogg	7.73	1st Learning	Non-target syllable	2510.0	1.73
17	p01	25SG.ogg	7.73	1st Learning	Non-target syllable	5079.0	1.73
18	p01	12SG.ogg	7.73	1st Learning	Non-target syllable	2622.0	1.73
19	p01	23PA.ogg	7.73	1st Learning	Target syllable	2161.0	1.73
20	p01	32PA.ogg	7.73	1st Learning	Target syllable	2141.0	1.73
21	p01	19PA.ogg	7.73	1st Learning	Target syllable	2311.0	1.73
22	p01	13PA.ogg	7.73	1st Learning	Target syllable	3049.0	1.73
23	p01	04PA.ogg	7.73	1st Learning	Target syllable	2175.0	1.73
24	p01	25PA.ogg	7.73	1st Learning	Target syllable	1750.0	1.73
25	p01	20PA.ogg	7.73	1st Learning	Target syllable	2158.0	1.73
26	p01	14PA.ogg	7.73	1st Learning	Target syllable	3274.0	1.73
27	p01	11PA.ogg	7.73	1st Learning	Target syllable	2434.0	1.73
28	p01	30PA.ogg	7.73	1st Learning	Target syllable	1859.0	1.73
29	p01	15PA.ogg	7.73	1st Learning	Target syllable	2351.0	1.73
30	p01	16PA.ogg	7.73	1st Learning	Target syllable	1853.0	1.73
⋮	⋮	⋮	⋮	⋮	⋮	⋮	⋮

• Centering age at six years yields an interpretable GM
• Factor levels can also be set when contrasts are defined (see below).
• BoxCox check showed that reaction time rt [ms] should be transformed to speed [1/s] = [Hz]
• Indicator variables for Target and Block generated in R.

4 LMM analysis

4.1 Contrasts

Code

contrasts = merge(
      Dict(:Target => EffectsCoding()),
      Dict(:Block => SeqDiffCoding()),
      Dict(nm => Grouping() for nm in (:Subj, :Item))
   );

4.2 Varying only intercepts LMM m_voi

Code

f_voi1 = @formula(lrt ~  1 + Block * Target * a1 + (1 | Subj) + (1 | Item));
m_voi1 = fit(MixedModel, f_voi1, dat; contrasts)

Minimizing 48    Time: 0:00:01 (24.56 ms/it)

	Est.	SE	z	p	σ_Item	σ_Subj
(Intercept)	7.7048	0.0358	215.08	<1e-99	0.0255	0.2508
Block: 2nd Learning	-0.0810	0.0102	-7.95	<1e-14
Block: Disruption	0.0518	0.0124	4.18	<1e-04
Block: Recovery	-0.0509	0.0123	-4.14	<1e-04
Target: Target syllable	-0.0121	0.0040	-3.05	0.0023
a1	-0.1054	0.0314	-3.36	0.0008
Block: 2nd Learning & Target: Target syllable	-0.0003	0.0102	-0.03	0.9746
Block: Disruption & Target: Target syllable	-0.0383	0.0122	-3.15	0.0016
Block: Recovery & Target: Target syllable	0.0102	0.0121	0.85	0.3959
Block: 2nd Learning & a1	0.0237	0.0089	2.65	0.0081
Block: Disruption & a1	0.0211	0.0106	1.99	0.0470
Block: Recovery & a1	-0.0127	0.0106	-1.20	0.2305
Target: Target syllable & a1	-0.0040	0.0035	-1.15	0.2494
Block: 2nd Learning & Target: Target syllable & a1	-0.0165	0.0089	-1.85	0.0641
Block: Disruption & Target: Target syllable & a1	0.0234	0.0106	2.21	0.0274
Block: Recovery & Target: Target syllable & a1	0.0005	0.0105	0.04	0.9656
Residual	0.2919

4.3 Extract indicator variables

Code

X = modelmatrix(m_voi1)
dat.trng = X[:,2];
dat.drpt = X[:,3];
dat.rcvr = X[:,4];
dat.trgt = X[:,5];
describe(dat)

12 rows × 7 columns

	variable	mean	min	median	max	nmissing	eltype
	Symbol	Union…	Any	Union…	Any	Int64	Type
1	Subj		p01		p85	0	Union{Missing, String}
2	Item		01PA.ogg		32SG.ogg	0	Union{Missing, String}
3	age	6.23114	3.7	6.31	8.04	0	Union{Missing, Float64}
4	Block		1st Learning		Recovery	0	Union{Missing, CategoricalValue{String, UInt32}}
5	Target		Non-target syllable		Target syllable	0	Union{Missing, CategoricalValue{String, UInt32}}
6	rt	2343.97	228.0	2105.0	8389.0	0	Union{Missing, Float64}
7	a1	0.231143	-2.3	0.31	2.04	0	Float64
8	lrt	7.67966	5.42935	7.65207	9.03468	0	Float64
9	trng	-0.0140243	-0.75	0.25	0.25	0	Float64
10	drpt	-0.0518619	-0.5	-0.5	0.5	0	Float64
11	rcvr	0.0530206	-0.25	-0.25	0.75	0	Float64
12	trgt	0.00175803	-1.0	1.0	1.0	0	Float64

Switch to indicator variables and refit the model.

Code

f_voi2 = @formula(lrt ~  1 + (trng+drpt+rcvr) * trgt * a1 + (1 | Subj) + (1 | Item));
m_voi2 = fit(MixedModel, f_voi2, dat; contrasts)

	Est.	SE	z	p	σ_Item	σ_Subj
(Intercept)	7.7048	0.0358	215.08	<1e-99	0.0255	0.2508
trng	-0.0810	0.0102	-7.95	<1e-14
drpt	0.0518	0.0124	4.18	<1e-04
rcvr	-0.0509	0.0123	-4.14	<1e-04
trgt	-0.0121	0.0040	-3.05	0.0023
a1	-0.1054	0.0314	-3.36	0.0008
trng & trgt	-0.0003	0.0102	-0.03	0.9746
drpt & trgt	-0.0383	0.0122	-3.15	0.0016
rcvr & trgt	0.0102	0.0121	0.85	0.3959
trng & a1	0.0237	0.0089	2.65	0.0081
drpt & a1	0.0211	0.0106	1.99	0.0470
rcvr & a1	-0.0127	0.0106	-1.20	0.2305
trgt & a1	-0.0040	0.0035	-1.15	0.2494
trng & trgt & a1	-0.0165	0.0089	-1.85	0.0641
drpt & trgt & a1	0.0234	0.0106	2.21	0.0274
rcvr & trgt & a1	0.0005	0.0105	0.04	0.9656
Residual	0.2919

They are equivalent.

4.4 A zero-correlation parameter LMM m_zcp

Code

f_zcp1 = @formula(lrt ~  1 +  Block * Target * a1 + zerocorr(1 + Block * Target | Subj) + (1 + a1 | Item));
m_zcp1 = fit(MixedModel, f_zcp1, dat; contrasts);

show(issingular(m_zcp1))
VarCorr(m_zcp1)

Minimizing 375   Time: 0:00:00 ( 1.46 ms/it)

false

	Column	Variance	Std.Dev	Corr.
Subj	(Intercept)	0.0668442	0.2585425
	Block: 2nd Learning	0.0305112	0.1746747	.
	Block: Disruption	0.0165614	0.1286909	.	.
	Block: Recovery	0.0135351	0.1163403	.	.	.
	Target: Target syllable	0.0002921	0.0170918	.	.	.	.
	Block: 2nd Learning & Target: Target syllable	0.0003617	0.0190193	.	.	.	.	.
	Block: Disruption & Target: Target syllable	0.0023922	0.0489103	.	.	.	.	.	.
	Block: Recovery & Target: Target syllable	0.0006247	0.0249938	.	.	.	.	.	.	.
Item	(Intercept)	0.0009965	0.0315682
	a1	0.0005789	0.0240613	-0.77
Residual		0.0724795	0.2692202

Again, check the equivalence.

Code

f_zcp2 = @formula(lrt ~  1 + (trng+drpt+rcvr) * trgt * a1  +
                zerocorr(1 + (trng+drpt+rcvr) * trgt | Subj) + (1 + a1 | Item));
m_zcp2 = fit(MixedModel, f_zcp2, dat; contrasts);

show(issingular(m_zcp2))
VarCorr(m_zcp2)

Minimizing 375   Time: 0:00:00 ( 0.78 ms/it)

false

	Column	Variance	Std.Dev	Corr.
Subj	(Intercept)	0.0668442	0.2585425
	trng	0.0305112	0.1746747	.
	drpt	0.0165614	0.1286909	.	.
	rcvr	0.0135351	0.1163403	.	.	.
	trgt	0.0002921	0.0170918	.	.	.	.
	trng & trgt	0.0003617	0.0190193	.	.	.	.	.
	drpt & trgt	0.0023922	0.0489103	.	.	.	.	.	.
	rcvr & trgt	0.0006247	0.0249938	.	.	.	.	.	.	.
Item	(Intercept)	0.0009965	0.0315682
	a1	0.0005789	0.0240613	-0.77
Residual		0.0724795	0.2692202

4.5 A complex parameter LMM m_cpx

Code

m_cpx = let
    form = @formula(lrt ~  1 + trgt * (trng+drpt+rcvr) * a1  +
                          (1 + trgt * (trng+drpt+rcvr) | Subj) + (1 + a1 | Item));
    fit(MixedModel, form, dat; contrasts);
  end;

show(issingular(m_cpx))    # not ok
show(m_cpx.PCA[:Subj])     # not ok
show(MixedModels.likelihoodratiotest(m_zcp2, m_cpx))
VarCorr(m_cpx)

Minimizing 1640      Time: 0:00:01 ( 0.64 ms/it)
  objective:  1954.9673786896176

true

Principal components based on correlation matrix

 (Intercept)   1.0     .      .      .      .      .      .     .
 trgt         -0.34   1.0     .      .      .      .      .     .
 trng          0.23   0.13   1.0     .      .      .      .     .
 drpt         -0.01  -0.08   0.3    1.0     .      .      .     .
 rcvr         -0.35   0.01  -0.13  -0.45   1.0     .      .     .
 trgt & trng   0.06  -0.52  -0.49   0.01   0.02   1.0     .     .
 trgt & drpt  -0.16   0.06   0.44   0.42  -0.38  -0.64   1.0    .
 trgt & rcvr   0.54   0.26  -0.36  -0.16  -0.09   0.13  -0.6   1.0

Normalized cumulative variances:

[0.33, 0.5515, 0.7424, 0.8685, 0.9466, 0.9934, 1.0, 1.0]

Component loadings
 

               PC1    PC2    PC3    PC4    PC5    PC6    PC7    PC8
 (Intercept)  -0.12   0.6    0.28  -0.44   0.08   0.11   0.34  -0.47
 trgt          0.14  -0.33   0.59   0.39  -0.21  -0.29   0.39  -0.29
 trng          0.43   0.13   0.14  -0.48  -0.51  -0.39  -0.07   0.36
 drpt          0.31   0.36  -0.2    0.43  -0.57   0.43  -0.06  -0.19
 rcvr         -0.21  -0.54  -0.08  -0.41  -0.38   0.51   0.3   -0.02
 trgt & trng  -0.43   0.2   -0.45   0.18  -0.24  -0.39   0.55   0.2
 trgt & drpt   0.57   0.03  -0.09   0.04   0.4    0.25   0.57   0.34
 trgt & rcvr  -0.37   0.24   0.55   0.18  -0.08   0.31  -0.02   0.6

Model Formulae
1: lrt ~ 1 + trng + drpt + rcvr + trgt + a1 + trng & trgt + drpt & trgt + rcvr & trgt + trng & a1 + drpt & a1 + rcvr & a1 + trgt & a1 + trng & trgt & a1 + drpt & trgt & a1 + rcvr & trgt & a1 + MixedModels.ZeroCorr((1 + trng + drpt + rcvr + trgt + trng & trgt + drpt & trgt + rcvr & trgt | Subj)) + (1 + a1 | Item)
2: lrt ~ 1 + trgt + trng + drpt + rcvr + a1 + trgt & trng + trgt & drpt + trgt & rcvr + trgt & a1 + trng & a1 + drpt & a1 + rcvr & a1 + trgt & trng & a1 + trgt & drpt & a1 + trgt & rcvr & a1 + (1 + trgt + trng + drpt + rcvr + trgt & trng + trgt & drpt + trgt & rcvr | Subj) + (1 + a1 | Item)
──────────────────────────────────────────────────
     model-dof  -2 logLik       χ²  χ²-dof  P(>χ²)
──────────────────────────────────────────────────
[1]         28  2007.1480                         
[2]         56  1954.9674  52.1807      28  0.0037
──────────────────────────────────────────────────

	Column	Variance	Std.Dev	Corr.
Subj	(Intercept)	0.0676775	0.2601489
	trgt	0.0003092	0.0175844	-0.34
	trng	0.0288450	0.1698380	+0.23	+0.13
	drpt	0.0196619	0.1402209	-0.01	-0.08	+0.30
	rcvr	0.0180750	0.1344433	-0.35	+0.01	-0.13	-0.45
	trgt & trng	0.0015803	0.0397531	+0.06	-0.52	-0.49	+0.01	+0.02
	trgt & drpt	0.0065055	0.0806569	-0.16	+0.06	+0.44	+0.42	-0.38	-0.64
	trgt & rcvr	0.0030635	0.0553486	+0.54	+0.26	-0.36	-0.16	-0.09	+0.13	-0.60
Item	(Intercept)	0.0008437	0.0290467
	a1	0.0005576	0.0236127	-0.87
Residual		0.0720702	0.2684589

The deviance improves, but we end up with an overparameterized LMM.

4.6 A parsimonious parameter LMM m_prm

We remove one of the VC for trgt * trng contrast interaction, that is one of the three interaction terms.

Code

m_prm1 = let
    form = @formula(lrt ~  1 + trgt * (trng+drpt+rcvr) * a1  +
                          (1 + trgt *(drpt + rcvr) + trng  | Subj) +
                          (1 + a1 | Item));
    fit(MixedModel, form, dat; contrasts);
  end;

show(issingular(m_prm1))  # ok
show(m_prm1.PCA[:Subj])   # not ok
show(MixedModels.likelihoodratiotest(m_zcp2, m_prm1, m_cpx))
VarCorr(m_prm1)

Minimizing 2205      Time: 0:00:01 ( 0.63 ms/it)

true
Principal components based on correlation matrix
 (Intercept)   1.0     .      .      .      .      .      .
 trgt         -0.33   1.0     .      .      .      .      .
 drpt         -0.01  -0.11   1.0     .      .      .      .
 rcvr         -0.34   0.01  -0.45   1.0     .      .      .
 trgt & drpt  -0.18  -0.09   0.5   -0.43   1.0     .      .
 trgt & rcvr   0.54   0.28  -0.16  -0.09  -0.67   1.0     .
 trng          0.21   0.08   0.3   -0.12   0.4   -0.38   1.0

Normalized cumulative variances:
[0.3473, 0.5941, 0.7601, 0.8926, 0.9713, 1.0, 1.0]

Component loadings
                PC1    PC2    PC3    PC4    PC5    PC6    PC7
 (Intercept)  -0.1    0.68  -0.24   0.32  -0.11   0.3   -0.52
 trgt         -0.12  -0.16   0.87   0.19  -0.1    0.13  -0.37
 drpt          0.44   0.23   0.17  -0.28   0.78   0.02  -0.19
 rcvr         -0.29  -0.5   -0.26   0.34   0.46   0.51  -0.07
 trgt & drpt   0.58  -0.05   0.05  -0.11  -0.33   0.7    0.22
 trgt & rcvr  -0.46   0.45   0.28  -0.07   0.19   0.29   0.62
 trng          0.38   0.1    0.08   0.8    0.12  -0.24   0.34

Model Formulae
1: lrt ~ 1 + trng + drpt + rcvr + trgt + a1 + trng & trgt + drpt & trgt + rcvr & trgt + trng & a1 + drpt & a1 + rcvr & a1 + trgt & a1 + trng & trgt & a1 + drpt & trgt & a1 + rcvr & trgt & a1 + MixedModels.ZeroCorr((1 + trng + drpt + rcvr + trgt + trng & trgt + drpt & trgt + rcvr & trgt | Subj)) + (1 + a1 | Item)
2: lrt ~ 1 + trgt + trng + drpt + rcvr + a1 + trgt & trng + trgt & drpt + trgt & rcvr + trgt & a1 + trng & a1 + drpt & a1 + rcvr & a1 + trgt & trng & a1 + trgt & drpt & a1 + trgt & rcvr & a1 + (1 + trgt + drpt + rcvr + trgt & drpt + trgt & rcvr + trng | Subj) + (1 + a1 | Item)
3: lrt ~ 1 + trgt + trng + drpt + rcvr + a1 + trgt & trng + trgt & drpt + trgt & rcvr + trgt & a1 + trng & a1 + drpt & a1 + rcvr & a1 + trgt & trng & a1 + trgt & drpt & a1 + trgt & rcvr & a1 + (1 + trgt + trng + drpt + rcvr + trgt & trng + trgt & drpt + trgt & rcvr | Subj) + (1 + a1 | Item)
──────────────────────────────────────────────────
     model-dof  -2 logLik       χ²  χ²-dof  P(>χ²)
──────────────────────────────────────────────────
[1]         28  2007.1480                         
[2]         48  1962.8756  44.2725      20  0.0014
[3]         56  1954.9674   7.9082       8  0.4425
──────────────────────────────────────────────────

	Column	Variance	Std.Dev	Corr.
Subj	(Intercept)	0.0678202	0.2604231
	trgt	0.0003055	0.0174778	-0.33
	drpt	0.0195723	0.1399012	-0.01	-0.11
	rcvr	0.0180836	0.1344754	-0.34	+0.01	-0.45
	trgt & drpt	0.0049548	0.0703907	-0.18	-0.09	+0.50	-0.43
	trgt & rcvr	0.0029429	0.0542485	+0.54	+0.28	-0.16	-0.09	-0.67
	trng	0.0297813	0.1725726	+0.21	+0.08	+0.30	-0.12	+0.40	-0.38
Item	(Intercept)	0.0008334	0.0288684
	a1	0.0005466	0.0233798	-0.87
Residual		0.0722829	0.2688548

We don’t lose goodness of fit, but are still overparameterized. We remove CPs for Target.

Code

m_prm2 = let
    form = @formula(lrt ~  1 + trgt * (trng+drpt+rcvr) * a1  +
                          (1 + drpt + rcvr + trng + drpt&trgt + rcvr&trgt  | Subj) +
                  zerocorr(0 + trgt  | Subj) +
                          (1 + a1 | Item));
    fit(MixedModel, form, dat; contrasts);
  end;

show(issingular(m_prm2))  #  ok
show(m_prm2.PCA[:Subj])   #  ok
show(MixedModels.likelihoodratiotest(m_zcp2, m_prm2,  m_cpx))
VarCorr(m_prm2)

Minimizing 1064      Time: 0:00:00 ( 0.58 ms/it)

false
Principal components based on correlation matrix
 (Intercept)   1.0     .      .      .      .      .     .
 drpt         -0.01   1.0     .      .      .      .     .
 rcvr         -0.34  -0.45   1.0     .      .      .     .
 trng          0.21   0.3   -0.12   1.0     .      .     .
 drpt & trgt  -0.12   0.53  -0.44   0.42   1.0     .     .
 rcvr & trgt   0.5   -0.2   -0.09  -0.44  -0.67   1.0    .
 trgt          0.0    0.0    0.0    0.0    0.0    0.0   1.0

Normalized cumulative variances:
[0.3523, 0.5915, 0.7343, 0.869, 0.9457, 0.9836, 1.0]

Component loadings
                PC1    PC2   PC3    PC4    PC5    PC6    PC7
 (Intercept)  -0.08   0.65  0.0   -0.44   0.14  -0.52  -0.3
 drpt          0.45   0.22  0.0    0.34  -0.76  -0.17  -0.18
 rcvr         -0.3   -0.54  0.0   -0.36  -0.41  -0.54   0.17
 trng          0.41   0.07  0.0   -0.73  -0.22   0.44   0.23
 drpt & trgt   0.58  -0.03  0.0    0.11   0.35  -0.46   0.56
 rcvr & trgt  -0.46   0.48  0.0    0.14  -0.25   0.09   0.69
 trgt          0.0    0.0   1.0    0.0    0.0    0.0    0.0

Model Formulae
1: lrt ~ 1 + trng + drpt + rcvr + trgt + a1 + trng & trgt + drpt & trgt + rcvr & trgt + trng & a1 + drpt & a1 + rcvr & a1 + trgt & a1 + trng & trgt & a1 + drpt & trgt & a1 + rcvr & trgt & a1 + MixedModels.ZeroCorr((1 + trng + drpt + rcvr + trgt + trng & trgt + drpt & trgt + rcvr & trgt | Subj)) + (1 + a1 | Item)
2: lrt ~ 1 + trgt + trng + drpt + rcvr + a1 + trgt & trng + trgt & drpt + trgt & rcvr + trgt & a1 + trng & a1 + drpt & a1 + rcvr & a1 + trgt & trng & a1 + trgt & drpt & a1 + trgt & rcvr & a1 + (1 + drpt + rcvr + trng + drpt & trgt + rcvr & trgt | Subj) + MixedModels.ZeroCorr((0 + trgt | Subj)) + (1 + a1 | Item)
3: lrt ~ 1 + trgt + trng + drpt + rcvr + a1 + trgt & trng + trgt & drpt + trgt & rcvr + trgt & a1 + trng & a1 + drpt & a1 + rcvr & a1 + trgt & trng & a1 + trgt & drpt & a1 + trgt & rcvr & a1 + (1 + trgt + trng + drpt + rcvr + trgt & trng + trgt & drpt + trgt & rcvr | Subj) + (1 + a1 | Item)
──────────────────────────────────────────────────
     model-dof  -2 logLik       χ²  χ²-dof  P(>χ²)
──────────────────────────────────────────────────
[1]         28  2007.1480                         
[2]         42  1969.2274  37.9206      14  0.0005
[3]         56  1954.9674  14.2600      14  0.4305
──────────────────────────────────────────────────

	Column	Variance	Std.Dev	Corr.
Subj	(Intercept)	0.0678705	0.2605197
	drpt	0.0195960	0.1399859	-0.01
	rcvr	0.0181233	0.1346227	-0.34	-0.45
	trng	0.0298994	0.1729145	+0.21	+0.30	-0.12
	drpt & trgt	0.0048302	0.0694998	-0.12	+0.53	-0.44	+0.42
	rcvr & trgt	0.0024143	0.0491356	+0.50	-0.20	-0.09	-0.44	-0.67
	trgt	0.0003211	0.0179201	.	.	.	.	.	.
Item	(Intercept)	0.0008544	0.0292300
	a1	0.0005640	0.0237488	-0.83
Residual		0.0723240	0.2689312

Check the Item-related CP. It is very large, might be spurious.

Code

m_prm3 = let
    form = @formula(lrt ~  1 + trgt * (trng+drpt+rcvr) * a1  +
                          (1 + drpt + rcvr + trng + drpt&trgt + rcvr&trgt  | Subj) +
                  zerocorr(0 + trgt  | Subj) +
                  zerocorr(1 + a1 | Item));
    fit(MixedModel, form, dat; contrasts);
  end;

show(issingular(m_prm3))  # ok
show(m_prm3.PCA[:Subj])   # ok
show(MixedModels.likelihoodratiotest(m_zcp2, m_prm3, m_prm2, m_cpx))
VarCorr(m_prm3)

Minimizing 1133      Time: 0:00:00 ( 0.55 ms/it)

false
Principal components based on correlation matrix
 (Intercept)   1.0     .      .      .      .      .     .
 drpt         -0.01   1.0     .      .      .      .     .
 rcvr         -0.34  -0.45   1.0     .      .      .     .
 trng          0.21   0.31  -0.11   1.0     .      .     .
 drpt & trgt  -0.13   0.53  -0.44   0.43   1.0     .     .
 rcvr & trgt   0.56  -0.17  -0.14  -0.45  -0.63   1.0    .
 trgt          0.0    0.0    0.0    0.0    0.0    0.0   1.0

Normalized cumulative variances:
[0.348, 0.5965, 0.7394, 0.8752, 0.9496, 0.9853, 1.0]

Component loadings
                PC1    PC2   PC3    PC4    PC5    PC6    PC7
 (Intercept)  -0.1    0.63  0.0   -0.46   0.1   -0.47  -0.39
 drpt          0.45   0.24  0.0    0.33  -0.77  -0.11  -0.19
 rcvr         -0.28  -0.54  0.0   -0.35  -0.44  -0.54   0.13
 trng          0.42   0.06  0.0   -0.73  -0.18   0.42   0.3
 drpt & trgt   0.58   0.0   0.0    0.12   0.35  -0.55   0.47
 rcvr & trgt  -0.45   0.49  0.0    0.13  -0.22   0.01   0.7
 trgt          0.0    0.0   1.0    0.0    0.0    0.0    0.0

Model Formulae
1: lrt ~ 1 + trng + drpt + rcvr + trgt + a1 + trng & trgt + drpt & trgt + rcvr & trgt + trng & a1 + drpt & a1 + rcvr & a1 + trgt & a1 + trng & trgt & a1 + drpt & trgt & a1 + rcvr & trgt & a1 + MixedModels.ZeroCorr((1 + trng + drpt + rcvr + trgt + trng & trgt + drpt & trgt + rcvr & trgt | Subj)) + (1 + a1 | Item)
2: lrt ~ 1 + trgt + trng + drpt + rcvr + a1 + trgt & trng + trgt & drpt + trgt & rcvr + trgt & a1 + trng & a1 + drpt & a1 + rcvr & a1 + trgt & trng & a1 + trgt & drpt & a1 + trgt & rcvr & a1 + (1 + drpt + rcvr + trng + drpt & trgt + rcvr & trgt | Subj) + MixedModels.ZeroCorr((0 + trgt | Subj)) + MixedModels.ZeroCorr((1 + a1 | Item))
3: lrt ~ 1 + trgt + trng + drpt + rcvr + a1 + trgt & trng + trgt & drpt + trgt & rcvr + trgt & a1 + trng & a1 + drpt & a1 + rcvr & a1 + trgt & trng & a1 + trgt & drpt & a1 + trgt & rcvr & a1 + (1 + drpt + rcvr + trng + drpt & trgt + rcvr & trgt | Subj) + MixedModels.ZeroCorr((0 + trgt | Subj)) + (1 + a1 | Item)
4: lrt ~ 1 + trgt + trng + drpt + rcvr + a1 + trgt & trng + trgt & drpt + trgt & rcvr + trgt & a1 + trng & a1 + drpt & a1 + rcvr & a1 + trgt & trng & a1 + trgt & drpt & a1 + trgt & rcvr & a1 + (1 + trgt + trng + drpt + rcvr + trgt & trng + trgt & drpt + trgt & rcvr | Subj) + (1 + a1 | Item)
──────────────────────────────────────────────────
     model-dof  -2 logLik       χ²  χ²-dof  P(>χ²)
──────────────────────────────────────────────────
[1]         28  2007.1480                         
[2]         41  1980.2755  26.8726      13  0.0130
[3]         42  1969.2274  11.0480       1  0.0009
[4]         56  1954.9674  14.2600      14  0.4305
──────────────────────────────────────────────────

	Column	Variance	Std.Dev	Corr.
Subj	(Intercept)	0.0681020	0.2609636
	drpt	0.0195454	0.1398050	-0.01
	rcvr	0.0181474	0.1347124	-0.34	-0.45
	trng	0.0299627	0.1730973	+0.21	+0.31	-0.11
	drpt & trgt	0.0043532	0.0659788	-0.13	+0.53	-0.44	+0.43
	rcvr & trgt	0.0019832	0.0445329	+0.56	-0.17	-0.14	-0.45	-0.63
	trgt	0.0003244	0.0180111	.	.	.	.	.	.
Item	(Intercept)	0.0007211	0.0268527
	a1	0.0004408	0.0209959	.
Residual		0.0724335	0.2691348

Perhaps not. We stay with LMM m_prm2

4.7 Compare models with goodness-of-fit statistics.

Code

let mods = [m_zcp2, m_prm2, m_prm1, m_cpx];
 DataFrame(;
    model=[:m_zcp2, :m_prm2, :m_prm1, :m_cpx],
    pars=dof.(mods),
    geomdof=round.(Int, (sum ∘ leverage).(mods)),
    AIC=round.(Int, aic.(mods)),
    AICc=round.(Int, aicc.(mods)),
    BIC=round.(Int, bic.(mods)),
  )
end

4 rows × 6 columns

	model	pars	geomdof	AIC	AICc	BIC
	Symbol	Int64	Int64	Int64	Int64	Int64
1	m_zcp2	28	293	2063	2063	2252
2	m_prm2	42	279	2053	2054	2336
3	m_prm1	48	274	2059	2060	2382
4	m_cpx	56	281	2067	2068	2444

LMM m_prm2 is a defensible solution according to \(\Delta\) AIC, \(\Delta\) BIC suggests we should not bother with CPs.

Code

coeftable(m_prm2)

	Coef.	Std. Error	z	Pr(>
(Intercept)	7.69869	0.0371901	207.01	<1e-99
trgt	-0.0118293	0.00451408	-2.62	0.0088
trng	-0.070356	0.0262888	-2.68	0.0074
drpt	0.0457635	0.0231034	1.98	0.0476
rcvr	-0.042934	0.0224154	-1.92	0.0554
a1	-0.108394	0.0326348	-3.32	0.0009
trgt & trng	0.00212125	0.00949196	0.22	0.8232
trgt & drpt	-0.0370438	0.0150283	-2.46	0.0137
trgt & rcvr	0.00972961	0.0132374	0.74	0.4623
trgt & a1	-0.00701276	0.00399025	-1.76	0.0788
trng & a1	0.0279914	0.0230904	1.21	0.2254
drpt & a1	0.0273081	0.0201717	1.35	0.1758
rcvr & a1	-0.0167551	0.0195894	-0.86	0.3924
trgt & trng & a1	-0.00608894	0.00831708	-0.73	0.4641
trgt & drpt & a1	0.029943	0.0131747	2.27	0.0230
trgt & rcvr & a1	-0.00509192	0.011614	-0.44	0.6611

4.8 Diagnostic plots

Various visualizations are used to check whether or not data are defensibly modeled with an LMM. They may lead to removal of outliers, transformations of the dependent variable, and deliver valuable heuristic information to be followed up with exploratory post-hoc analyses or ideally replication of new insights gained this way. In practice, it appears that only severe violations will stop people from reporting a model.

4.8.1 Residuals over fitted

Code

scatter(fitted(m_prm2), residuals(m_prm2))

Looks like we missed some fast response times.

4.8.2 Q-Q plot

Code

qqnorm(m_prm2; qqline=:none)

Hm?

4.8.3 Residual distributions: observed vs. theoretical

Curves for residulas based on observed and theoretical values should correspond.

Code

let
  n = nrow(dat)
  dat_rz = (;
    value=vcat(residuals(m_prm2) ./ std(residuals(m_prm2)), randn(n)),
    curve=repeat(["residual", "normal"]; inner=n),
  )
  draw(
    data(dat_rz) *
    mapping(:value; color=:curve) *
    density(; bandwidth=0.1);
  )
end

Figure 1: Kernel density plot of the standardized residuals for model m1 versus a standard normal

They are a bit too narrow.

4.9 Conditional means of random effects

4.9.1 Subject-related conditional means of random effects

Code

cm_Subj = first(ranefinfo(m_prm2))
caterpillar!(Figure(; resolution=(800, 1200)), cm_Subj; orderby=1)

Figure 2: Prediction intervals on subject random effects for model m_prm2

4.9.2 Borrowing-strength plots

Shrinkage refers to the adjustment of subject-level or item-level predictions by taking population estimates into account. The further a subject’s/item’s estimate is from the fixed effect or the more variable or less reliable the subject’s/item’s estimate, the more the prediction will be shrunk towards the population estimate. Alternative terms for shrinkage are “borrowing strength” (Tukey) and regularization. My favorite is actually Tukey’s because indeed we borrow strength from the population estimates to make predictions for individual subjects’ effects. The goal of this section to illustrate the results of borrowing strength.

Subject-related conditional means of random effects revealed information about individual differences beyond fixed effects. Would these results also be visible in unconditional means, that is when we compute GM and experimental effects within subjects (i.e., as fixed effects) without borrowing strength from the population estimates?

In the following plots, effect estimates based on alone on each subject’s data (i.e., no pooling of data, no borrowing of strength) are plotted in pink and the subjects’ conditional means shown in the caterpillar plots are plotted in blue. The arrows indicate how much a subject’s prediction is changed by borrowing strength from knowledge of the population estimates.

Code

shrinkageplot!(Figure(; resolution=(1000, 1200)), m_prm2)

Figure 3: Shrinkage plots of the subject random effects in model m_prm2

Outlier subject is gone.

Code

cm = raneftables(m_prm2);
sort(DataFrame(cm.Subj), ["(Intercept)", :trng, "drpt & trgt"])
sort(DataFrame(cm.Subj), :trng)

51 rows × 8 columns (omitted printing of 1 columns)

	Subj	(Intercept)	drpt	rcvr	trng	drpt & trgt	rcvr & trgt
	String	Float64	Float64	Float64	Float64	Float64	Float64
1	p14	-0.528323	0.159662	-0.118234	-0.385548	0.0520919	-0.00337878
2	p76	0.064757	-0.095427	0.0965818	-0.317894	-0.061881	0.0540763
3	p36	0.373846	-0.0446188	0.114011	-0.277302	-0.12993	0.109552
4	p02	-0.628903	-0.185906	0.239485	-0.24184	-0.0114483	-0.0500986
5	p39	0.292941	-0.123337	0.131867	-0.189144	-0.0662858	0.0443946
6	p74	-0.334411	-0.035486	0.058828	-0.177365	-0.00145978	-0.0244133
7	p67	0.0102658	0.00522903	-0.0258083	-0.171838	0.00596919	0.0146027
8	p59	0.376658	0.0232177	-0.253912	-0.137989	0.00689185	0.0472082
9	p06	-0.238811	0.140842	0.0320315	-0.128487	-0.0423011	0.0147919
10	p54	0.319235	-0.162324	0.0453531	-0.110458	-0.045112	0.0524742
11	p42	-0.105171	0.295914	-0.174955	-0.107408	0.0974469	-0.00806991
12	p23	-0.0173974	-0.0770885	0.064224	-0.105342	-0.0922702	0.0454142
13	p30	-0.34938	-0.153186	-0.0891651	-0.098345	0.00105422	-0.0225567
14	p45	-0.374579	-0.294247	0.154624	-0.089149	-0.00555676	-0.045438
15	p84	-0.432463	-0.215894	0.466653	-0.0825766	-0.170963	-0.013087
16	p83	-0.23897	-0.0537439	-0.162418	-0.0772424	0.0230343	-0.0184981
17	p72	0.177019	-0.0778771	-0.0337592	-0.0503027	-0.0629824	0.0679022
18	p55	-0.24107	-0.0677181	-0.0693261	-0.0486936	-0.0284252	-0.0116131
19	p27	-0.013248	-0.0111215	0.0152268	-0.0443757	-0.0441275	0.0253996
20	p12	-0.160657	-0.0495452	0.0628843	-0.0434829	-0.035482	0.00879487
21	p69	0.0806806	-0.0491015	-0.0696128	-0.0317218	0.044258	0.025753
22	p04	0.141053	-0.0468601	0.0393693	-0.0239661	-0.0124761	0.0138384
23	p77	0.150813	-0.101454	0.0587802	-0.0141855	-0.0636443	0.0327982
24	p79	-0.0709235	-0.0186446	0.0306235	-0.0134209	-0.022047	0.000434288
25	p53	-0.112958	0.0249321	0.0212246	-0.00152289	0.0380108	-0.0234996
26	p08	-0.160325	0.11263	0.0336374	0.00561131	0.0106376	-0.000935882
27	p25	0.278466	0.141361	-0.0476478	0.00933614	-0.0249672	0.0311902
28	p64	0.266069	0.00293777	0.034395	0.0113388	-0.0173697	0.0242711
29	p28	0.151958	-0.0047019	-0.159594	0.0136686	0.0485822	-0.00713976
30	p57	0.00152955	-0.0017727	-0.0199141	0.0196826	0.037071	-0.01744
⋮	⋮	⋮	⋮	⋮	⋮	⋮	⋮

4.9.3 Item-related conditional means of random effects

Code

cm_Item = last(ranefinfo(m_prm2))
caterpillar!(Figure(; resolution=(800, 1200)), cm_Item; orderby=1)

Figure 4: Prediction intervals on item random effects for model m_prm2

5 Appendix

Code

versioninfo()

Julia Version 1.8.1
Commit afb6c60d69a (2022-09-06 15:09 UTC)
Platform Info:
  OS: Linux (x86_64-linux-gnu)
  CPU: 8 × 11th Gen Intel(R) Core(TM) i5-1135G7 @ 2.40GHz
  WORD_SIZE: 64
  LIBM: libopenlibm
  LLVM: libLLVM-13.0.1 (ORCJIT, tigerlake)
  Threads: 8 on 8 virtual cores